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Abstract 

An isospin lattice-gas model, which is a spin-1 Ising model, is employed to 
investigate the liquid- gas phase transition in asymmetric nuclear matter. We 
consider nuclear matter as a lattice where each lattice site can be either empty 
or occupied by a proton or a neutron, with a nearest-neighbor interaction among 
the nucleons. With the Bragg- Williams mean field approximation, we calculate 
various thermodynamic properties of nuclear matter for different densities and 
different proton- neutron asymmetry parameter s. Our model exhibits liquid- 
gas phase transition below a critical temperature T c , and predicts a monotonic 
decreasing of T c as the magnitude of s is increased. The dependence of the 
nuclear matter isotherms on the asymmetry parameter s is discussed. 
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The study of the equation of state (EOS) of asymmetric nuclear matter has been, 
in the last few years, a subject of renewed interest particularly in connection with 
astrophysical problems!], 0, such as supernovae explosions and the evolution of neutron 
stars. The EOS for the latter samples a range of densities and isospin asymmetry which 
are different from those for supernovae. Thus these two physical systems provide a 
unique laboratory where the EOS of nuclear matter can be critically investigated and 
give the possibility of obtaining related but not identical information about the EOS. 

Several laboratories have studied the nuclear matter EOS by way of heavy-ion 
collisions, such as the recent 197 Au+ 197 Au reaction || investigated at GSI-Darmstadt. 
The proton-neutron ratio in 197 Au is approximately 2/3. Clearly, the nuclear matter 
formed via such reactions should be significantly asymmetric, lending further support 
to the need for carrying out theoretical studies of the EOS for asymmetric nuclear 
matter. 

In the past, asymmetric nuclear matter calculations have been done using standard 
many-body methods, such as the HF calculations using the Skyrme interactions 0, [5|] 
and the Gogny interaction ||. Brueckner-Hartree-Fock calculations for asymmetric 
nuclear matter, using realistic nucleon-nucleon (NN) interactions have also been per- 
formed 0, g, |. 

We would like to consider a different approach for the EOS of assymmetric nuclear 
matter. In statistical physics, phase transitions in extended systems are generally 
studied using lattice (Ising) models. For instance, the lattice gas model of Lee and 
Yang DTOf, where a gas with one type of atoms is mapped into an Ising model for spin-| 
particles, has successfully described the liquid-gas phase transitions for atomic systems. 
It should be of interest to explore the feasibilty of describing the nuclear-matter EOS 
using lattice-gas model. Indeed, several authors O, O O, [Tj|, [T5|] have employed 
lattice-gas models to study the liquid-gas phase transition in symmetric nuclear matter. 
In this paper we present a lattice-gas model for investigating the liquid-gas phase 
transitions in asymmetric nuclear matter. 
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We revisit the lattice hamiltonian for nuclear matter proposed by Kuo et.al.|TTl 115 



Hint ^ t Jij7~zi7~zj h ^ T 2 j 

<ij> i 



The summation index < ij > implies nearest neighbor interaction and h is some applied 
external field. are interaction strength parameters : J s for symmetric pairs (pp and 
nn) and Jd for asymmetric type pairs (pn pairs). The have the property 



for nearest neighbor distance = a 
Jij = \ oo for neighbor distance = (2) 
otherwise 



where a is the lattice spacing. The isospin r z is 



— 1 for neutron 
t z = ^ for vacancy (3) 

1 for proton 



If N ++ ,N and iV + „ represent respectively the nearest neighbor pairs of proton- 
proton, neutron-neutron and proton-neutron, the interaction hamiltonian may be writ- 
ten as, 

H int = ~Js(N ++ + iV__) + J d N + _ - h{N + - AL) (4) 
There is no interaction between vacancy and nucleon. 

3. Thermodynamics 

Let iV denote the total number of lattice sites, and iV + , iV_, N the number of proton, 
neutron and vacancy sites respectively. The relative emptiness (r) and proton-neutron 
asymmetry parameter (s) are defined as, 
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The nuclear-matter density p is proportional to (1 — r) in our model. 

For atomic systems, where the spin-| lattice gas models have been used with re- 



markable success in describing phase transitions, Huang |L6| suggests adding the ideal 
gas pressure to the lattice gas grand potential in order to describe the isotherms. The 
ideal gas pressure comes from the kinetic energy. This motivated us to include a 
kinetic-energy term in our hamiltonian itself, namely 

H gas = -J S (N ++ + iV__) + J d N + _ + Nk(1 - rf/ 3 - Nhs (6) 

where we have assumed the kinetic energy per particle to be proportional to p 2 / 3 as 
per the Fermi gas model, k is a constant which will be discussed later. From the 
above hamiltonian we can write down the partition function. For the evaluation of the 
partition function, we have employed the Bragg- Williams mean field approximation |nj, 
namely 

Nl N ++ N 2 iV__ N 2 iV 00 _ 
N 2 ~ AT 7 /2' N 2 ~ N 2 ~ Nj/2 v ; 

where 7 denotes the number of nearest neighbors to any given site, and is the 

total number of pairs. For three dimensional simple cubic lattice, 7 = 6. If one counts 
the number of bonds between a site and all its nearest neighbors, then for (say) all the 
proton sites one obtains 

1 N + = 2N ++ + iV + _ + iV +0 (8) 

and similarly for the neutron and vacant sites. 
Using the above we can eliminate iV + _ as, 

iV + _ = 7 (iV + + AL) - (N ++ + AL_) + N 00 - jN/2 (9) 

The hamiltonian of eq. (6) can then be rewritten as 

H gas (r, s, N) = -dNs 2 - C 2 N{1 - rf + Nk{1 - r) 5/3 - hNs (10) 



where 



(11) 
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Unlike the symmetric case, s is now a parameter and not a variable. The canonical 
partition function is therefore given by 

Q 9 as = J2 9(r,s,N)e^ H ^ (12) 



r=0 



where the multiplicity factor |L7| g(r, s, N) is given by 

9{r,S,N) = N \N + \NJ = (jVr)![f(l-r + s)]![f(l-r-s)]! (13) 
To study the system at varying densities, we need to work with the grand partition 

function with a chemical potential A 



Qg = E 9(r,s,N) 



(N++N^) J-f3H gas ) 



= E 9(r, s, N) Z ( N + +N ^ e (- 0H -) (14) 

r 

where z = e^ x is the fugacity. 

In the thermodynamic limit (N — > oo), the sum in eq. (14) can be replaced by 
its most dominant term [H] (assuming the dominant term to be non-degenerate). We 
have verified that this approximation is acceptable by numerically checking the relative 
size of the different leading order terms. Using Stirling's formula one obtains 

In Q G = pC x Ns* + (3C 2 N(1 - rf + (3hNs + (3XN(1 - r) - (3kN(1 - r) 5/3 
_ jV[ (l_-r+_s) ln(l _ r + s) + (1-r-s) ln(1 _ r _ s) 

+rlnr- (1 -r)ln2] (15) 
with the extremum condition 

+/5A^A = (16) 



In a previous work [15| it was shown that h = and s = results in a global maximum 
for this term. In the present case however we allow s to vary parametrically. Factoring 
out the common terms eq. (16) is written as 

k B T ln[- - 2 -] + 2C 2 {1 - r) - §«(1 - r) 2 / 3 + A = 0. (17) 

1 — r — s z 6 

A graphical analysis of the above equation reveals the existence of a condensed phase. 
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4. Phase transition and critical temperature 

Let us consider the A=0 case first. In this case we rewrite eq. (17) as 

x(r,T) = k B T-f(r-T)-g( r ) = (18) 

with 

1 4r 2 
/W " 2 ln[ (l-r) 2 - s2 ] 

g(r) = y(l-r) 2 / 3 -2C 2 (l-r) 

In the region r G [0, 1], k B T ■ f(r) is a monotonically increasing, unbounded function 
of r. It has one point of inflection where it intersects the r axis and this is at r = | 
when s = 0. In the same domain g{r) however, is a bounded function with a negative 
curvature and one maximum in the region r G [0, 1]. Thus there is a possibility of g(r) 
intersecting with k B T • /(r) more than once below a suitable temperature T c . k B T- f(r) 
intersects the r axis at r = i+V 1 ^ 1 s ^ which depends on s but is independent of 
the value of T. On the other hand g(r) intersects the r axis at two points: One is 
at r = 1 which is independent of the values of the parameters C*2, n, and s, and the 
second point of intersection depends on s and the ratio a = 5/t/6C2. For a suitable 
choice of the parameter a, the functions k B T ■ f(r) and g(r) intersect the r axis at a 
common point given by r = ^V 1 ^ 1 s \ ^ j g then determined by that value 
of temperature for a given s at which the curves k B T ■ f(r) and g(r) are tangential at 
their common zero point (point of intersection with the r axis). 

In Fig. 1 we display some typical behaviours of k B T ■ /(r) and g(r). The points 
of intersection of the two curves are shown for a family of k B T ■ f(r) and g(r) with 
k B T = 5, 10, 15 MeV. This is done for several values of the asymmetry parameter. 

For T < T c , the curves k B T ■ f(r) and g(r) have three intersection points, as 
denoted by A, B, and C for the k B T = 5 MeV case. It is readily checked that the 
middle intersection point, i.e. B, corresponds to a minimum of ln^G an d hence it is 
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not the physically interesting solution that we are looking for. The intersection points 
A and C are the physical solutions. 

At some critical temperature T c the curves k B T c ■ f(r) and g(r) are tangent to 
each other at their common zero point. Above this critical temperature T c , the two 
curves intersect only once at their common point of intersection. Below T c however, 
the slope of k B T • /(r) at its zero point is less than that of g(r) and two more points 
of intersection appear which represent the densities of the two different phases. Thus 
the critical vacancy ro and the critical temperature T c are determined by 

k B T c ■ f(r ) = g(r ) 

k B T c -f(r ) = g'(r ) (19) 

As the magnitude of the asymmetry parameter s increases, the curve T • /(r) becomes 
steeper and multiple intersection with g(r) becomes possible only at lower tempera- 
tures. For pure neutron matter (s = —1.0) the critical temperature from our model is 
equal to and there is only one phase as is expected. (Fig. 3) 

5. Equation of state and isotherms 

The equation of state is obtained as 

P(r,T) = ^lnQ G 

= C lS 2 + hs + C 2 (l - ff + A(l - f) - k{1 - f) 5/3 

-k B T[ {1 ~ f 2 + S) ln(l -f + s)+ (1 ~^~ S) ln(l -f-s) 
+flnf- (1 -f)ln2] (20) 

i - -(i-o pd 

Note that the first two terms contribute just an additive constant to the pressure for 
any given s. We calculate the p-v isotherms using eqs.(18), (20) and (21). Our results 
are shown in Figs. 4 to 6. 

For a given temperature T the isotherm is obtained as a parametric plot of specific 
volume versus pressure with A its generating parameter. For T > T c one gets a single 
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smooth curve by varying A. The phase boundary corresponds to A = 0. The isotherms 
to the left of this boundary are obtained with positive A, those to the right with negative 
A. The figure suggests the existence of three regions below T c , the dense (small r, large 
p) liquid-like phase, the rare (large r, sma 11 p) gas-like phase and the coexistence phase 
in between. 

For T below T c , it is of interest to note that no isotherms are obtained in the 
intermediate region, i.e. the coexistence region. The above is because when T < T c , 
eq. (18) has no solutions in the coexistence region. With the introduction of an 
infinitesimal A the system chooses one of the two values of f (or v = nz?y ) admissible 
for the given T < T c , depending on the sign of A. This phenomena is reminiscent of 
the spontaneous symmetry breaking in ferromagnetism. 

The isotherms for different asymmetry parametered systems are shown in Fig. 4, 
5 and 6. These graphs are obtained by a very densely chosen set of points where eq. 
(18) is solved and the corresponding pressure and volume determined. To avoid falling 
into the trap of seeing a smooth isotherm when it is not, we intentionally chose not to 
use spline, or any higher polynomial fit. A very densely packed set of points with short 
straight line connects generated the isotherms. The absence of jaggedness that can be 
seen is coming from the physics and calculation itself, and not due to an artifact of 
any particular curve fitting technique. 

6. Analysis and results 

To have an attractive nearest-neighbor interaction, J s > and Jd < 0. In this case Ci 
is positive. The magnitude of J s and Jd is comparable to the average potential energy 
in nuclear matter, about —40 MeV. Hence C% ~ 125 MeV, recalling that Ci has been 
defined in eq. (11). The parameter k may be estimated from the average kinetic energy 
given by the Fermi gas model. 

Near the critical point, there is a subtle dependence of the solutions of eq. (18) on 
the ratio a = J=jn. For a choice of a = (2/3) 1 / 3 , the three intersection points below T c 
(A, B and C of Fig. 1), all merge together at the critical point as in the s = case|15 . 
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Then the phase diagram near T c has the familiar smooth shape. For a different choice of 
a, the merging would generally take place in two steps, first involving two intersection 
points and then the third. This will lead to a phase diagram with a cusp shape near the 
critical point, which seems rather unconventional. It would be interesting to investigate 
these unconventional phase boundaries further for different values of a. For the present 
study however, we have chosen a so as to get a smooth phase boundary. This choice 
of a is one for which the curves k^T ■ f(r) and g(r) intersect the r axis at a common 
point. With this a and C*2=125, we have k=131.037, which has been used in obtaining 
the results presented in our Figs. 1 to 6. 

For fixed T, the solutions of eq. (17) determine the r values where the grand 
partition function is maximum. Very similar to what was observed in the s = case, 
this maximum term is an overwhelmingly dominant one. (Numerical simulation for 
finite lattices has shown that the magnitude of the maximum term is typically a factor 
of ~ 10 80 greater than the other terms.) Hence the r values given by eq. (18) are just 
the average values (stricly speaking the modal value) of r, denoted by f, for the system 
at temperature T, asymmetry s and A=0 for which lnQc is maximum. The f values 
for A 7^ are given by the solution of eq. (17). 

We find support of our earlier observation |~5] in our present calculation as well, 
that a term of Fermi gas type kinetic energy is needed for appropriate phase structure 
to appear. If there is no kinetic energy term, k = and hence a = 0. Then the 
requirement that at T < T c , g{r) intersect T ■ f(r) more than once cannot be met, and 
our model would have no phase transitions. 

Conclusion 

Our simple model exhibits the existence of a liquid gas phase transition in asymmetric 
nuclear matter. The p-v isotherms obtained by our model look surprisingly similar 
to those given by the van der Waals theory, except for one crucial difference that for 
T < T c our isotherms do not have the metastable states in the coexistence region. 
Hence with our model one does not need to determine the phase boundary by way of 
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a Maxwell construction. 

The liquid-gas critical temperature obtained is dependent on the asymmetry pa- 
rameter s, T c decreases with increasing neutron density as shown in Fig. 3. The value 
is around T c = 17 — 20 for typical values of s which is fairly close to the results given 
by earlier calculations 0, [18|, [L9|, [y|. The existence of a liquid-gas phase transition 
together with the determination of its critical temperature on the basis of a simple 
model that assumes only a phenomenological two-body, nearest-neighbor interaction 
is quite remarkable. 



Ray, Shamanna and Kuo / Liquid-gas phase transition in Nuclear matter 11 

References 

[1] E. Baron, J. Cooperstein and S. Kahana, Phys. Rev. Lett. 55 (1985) 126. 
[2] M. Brack, C. Guet and H. B. Hakansson, Phys. Rep. 123 (1985) 277. 
[3] J. Pochodzalla et al, Phys. Rev. Lett. 75 (1995) 1040. 
[4] H. Jaqaman, A. Z. Mekjian and L. Zamick, Phys. Rev. C27 (1983) 2782. 
[5] R. K. Su, S. D. Yang, G. L. Li and T. T. S. Kuo, Mod. Phys. Lett. Al (1986) 71. 
[6] H. Q. Song, G. D. Zheng and R. K. Su, J. Phys. G16 (1990) 1861. 
[7] I. Bombaci and U. Lombardo, Phys. Rev. C44 (1991) 1892. 
[8] I. Bombaci, T. T. S. Kuo and U. Lombardo, Phys. Rep. 242 (1994) 165-180. 
[9] H. Q. Song, Z. X. Wang and T.T.S. Kuo, Phys. Rev. C46 (1992) 1788. 
[10] T.D. Lee and C.N. Yang, Phys. Rev. 87 (1952) 410. 

[11] R. K. Su and T. T. S. Kuo, Isospin Lattice Gas Model for Nuclear Matter, Stony 
Brook preprint (1990), unpublished. 

[12] J. Shamanna, S. Ray, T.T.S. Kuo, R.K. Su, Bull. Am. Phys. Soc. 39 (2) 1138 
(1994). 

[13] J. Pan and S. Das Gupta, Phys Lett B344 (1995) 29. 
[14] J. Pan and S. Das Gupta, Phys Rev C51 (1995) 1384. 

[15] T.T.S. Kuo, S. Ray, J. Shamanna, R.K. Su, Int. J. Mod. Phys. E5 No. 2 (1996) 
303. 

[16] K. Huang, Statistical Mechanics (Wiley, New York, 1963). 

[17] William Feller, Probability Theory (Wiley InterScience, 1976). 

[18] R. K. Su, S. D. Yang and T. T. S. Kuo, Phys. Rev C35 (1987) 1539. 



Ray, Shamanna and Kuo / Liquid-gas phase transition in Nuclear matter 



12 



[19] M. Baldo, I. Bombaci, L. S. Ferreira, G. Giansiracusa and U. Lombardo, Phys. 
Lett. B215 (1988) 19. 



Ray, Shamanna and Kuo / Liquid-gas phase transition in Nuclear matter 13 

FIGURE CAPTIONS 

Fig.l Graphical solution of eq.(18) for s=-0.2. 

Fig. 2 Graphical solution of eq.(18) for s=-0.4. 

Fig. 3 Critical Temperature versus asymmetry parameter. 

Fig. 4 Nuclear matter p-v isotherms for s=-0.2. 

Fig. 5 Nuclear matter p-v isotherms for s=-0.4. 

Fig. 6 Nuclear matter p-v isotherms for s=-0.6. 
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